import matplotlib.pyplot as plt
import numpy as np

from scipy.io import loadmat

def load_image(path_img):
    im_loaded = loadmat('./'+ path_img +'.mat')
    return im_loaded['imm']

# First image: SPOT
path_img = 'SPOT/cam'
n_bands = 3

# Let's take a look at the first channel's size
img = load_image(path_img +'1')
img_shape = img.shape
print(img_shape)

(512, 512)

img = np.zeros((img.shape[0], img.shape[1], n_bands))
# Stacking up images into the numpy array:
for i in range(n_bands):
    img[:,:,i] = load_image(path_img + str(i+1))

import matplotlib.pyplot as plt
import matplotlib.gridspec as grid

fig,axes = plt.subplots(1,n_bands, figsize=(n_bands*10,12), sharex='all', sharey='all')
fig.subplots_adjust(wspace=0.1, hspace=0.15)
axes = axes.ravel()

for i in range(n_bands):
    # Use floats instead of 8-bit integers for a clearer view of the image
    axes[i].imshow(img[:,:,i] / 255, cmap='gray')
    axes[i].set_title('band '+str(i+1),fontsize=15)
    axes[i].axis('off')

# Reshape the image for practicality: 
print("Starting shape: ", img.shape)
img = img.reshape(-1, n_bands)
print("New shape: ", img.shape)
print(img)

Starting shape:  (512, 512, 3)
New shape:  (262144, 3)
[[29. 23. 29.]
 [30. 24. 29.]
 [32. 27. 30.]
 ...
 [39. 28. 35.]
 [35. 23. 37.]
 [26. 17. 37.]]

# Standardization:
means = np.mean(img, axis=0)
stds = np.std(img, axis=0)

img = (img - means) / stds

print(img)

[[-0.33621727  0.31463317  0.12187387]
 [-0.06011854  0.5357671   0.12187387]
 [ 0.49207891  1.19916889  0.26818369]
 ...
 [ 2.42476999  1.42030282  0.99973279]
 [ 1.32037509  0.31463317  1.29235244]
 [-1.16451345 -1.01217042  1.29235244]]

# Compte the covariance matrix:
cov = np.cov(img, rowvar=False)

# Compute eigenvalues and vectors with eig:
eig_vals, eig_vecs = np.linalg.eig(cov)

# Check that that it worked using the formula:
# A . v = eigenvalues . v

print("A . v\t\t\t\t\teigenvalues . v\t\t\t\tEqual?")
for i in range(n_bands):
  result1 = np.dot(cov, eig_vecs[:, i])
  result2 = eig_vals[i] * eig_vecs[:, i] 
  is_equal = np.allclose(result1, result2)

  print(f"{result1}\t{result2}\t{is_equal}")

A . v					eigenvalues . v				Equal?
[-1.40100567 -1.44503693 -1.0857528 ]	[-1.40100567 -1.44503693 -1.0857528 ]	True
[-0.05498694  0.05913769 -0.00775418]	[-0.05498694  0.05913769 -0.00775418]	True
[-0.25691484 -0.16637967  0.55294714]	[-0.25691484 -0.16637967  0.55294714]	True

# Lets sort the eigenvalues to be able to choose the components by importance:
idx = np.argsort(eig_vals)[::-1]
eig_vals = eig_vals[idx]
eig_vecs = eig_vecs[:, idx]

# Lets get the components in separate matrices, so we can display them easily:
img_principal = np.dot(img, eig_vecs)

plt.figure(figsize=(n_bands*10, 12))

# Visualize them with imshow: 
for i in range(n_bands):
    plt.subplot(1,n_bands,i+1)
    plt.imshow(img_principal[:,i].reshape(img_shape[0], img_shape[1]), cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

# Compute the percentage of inertia of each component:
inertia = eig_vals / np.sum(eig_vals)

print(inertia)

[0.76228952 0.21066955 0.02704093]

path_img = 'Landsat_Tarascon/landsattarasconC'
n_bands = 8

# Let's take a look at the first channel's size
img = load_image(path_img + '1')
img_shape = img.shape

img = np.zeros((img.shape[0], img.shape[1], n_bands))

# Stacking up images into the numpy array:
for i in range(n_bands):
    img[:, :, i] = load_image(path_img + str(i+1))

# Visualize them with imshow: 
n_rows = 2
plt.figure(figsize=(n_bands*10/n_rows, 20))

for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.imshow(img[:, :, i] / 255, cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

# Reshape the image for practicality: 
img = img.reshape(-1, n_bands)

# Standardization:
means = np.mean(img, axis=0)
stds = np.std(img, axis=0)

img = (img - means) / stds

# Compte the covariance matrix:
cov = np.cov(img, rowvar=False)

# Compute eigenvalues and vectors with eig:
eig_vals, eig_vecs = np.linalg.eig(cov)

# Lets sort the eigenvalues to be able to choose the components by importance:
idx = np.argsort(eig_vals)[::-1]	
eig_vals = eig_vals[idx]
eig_vecs = eig_vecs[:, idx]

# Lets get the components in separate matrices, so we can display them easily:
img_principal = np.dot(img, eig_vecs)

plt.figure(figsize=(n_bands*10/n_rows, 20))

for i in range(n_bands):
    plt.subplot(n_rows,n_bands//n_rows,i+1)
    plt.imshow(img_principal[:,i].reshape(img_shape[0], img_shape[1]), cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

# Computing the percentage of inertia of each component:
inertia = eig_vals / np.sum(eig_vals)

print('Intertia:')
print(inertia)

objective = 0.95
num_of_images = 0

while np.sum(inertia[:num_of_images]) < objective:
    num_of_images += 1

print("")
print(f'Number of images to reach {objective*100}% of inertia: {num_of_images}')

Intertia:
[0.68529996 0.16981093 0.09404482 0.02665519 0.01677483 0.00364708
 0.00249158 0.00127561]

Number of images to reach 95.0% of inertia: 4

# Capteur normal:
path_img = 'Landsat_Kedougou/landsatKedougouC'
n_bands = 6

# Let's take a look at the first channel's size
img = load_image(path_img + '1')
img_shape = img.shape

img = np.zeros((img_shape[0], img_shape[1], n_bands))

# Stacking up images into the numpy array:
for i in range(n_bands):
    img[:, :, i] = load_image(path_img + str(i+1))

# Visualize them with imshow: 
n_rows = 2
plt.figure(figsize=(n_bands*10/n_rows, 20))

for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.imshow(img[:, :, i] / 255, cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

# Reshape the image for practicality: 
img = img.reshape(-1, n_bands)

# Standardization:
means = np.mean(img, axis=0)
stds = np.std(img, axis=0)

img = (img - means) / stds

# Compte the covariance matrix:
cov = np.cov(img, rowvar=False)

# Compute eigenvalues and vectors with eig:
eig_vals, eig_vecs = np.linalg.eig(cov)

# Lets sort the eigenvalues to be able to choose the components by importance:
idx = np.argsort(eig_vals)[::-1]
eig_vals = eig_vals[idx]
eig_vecs = eig_vecs[:, idx]

# Lets get the components in separate matrices, so we can display them easily:
img_principal = np.dot(img, eig_vecs)

plt.figure(figsize=(n_bands*10/n_rows, 20))

for i in range(n_bands):
    plt.subplot(n_rows,n_bands//n_rows,i+1)
    plt.imshow(img_principal[:,i].reshape(img_shape[0], img_shape[1]), cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

# Computing the percentage of inertia of each component:
inertia = eig_vals / np.sum(eig_vals)

print('Intertia:')
print(inertia)

objective = 0.95
num_of_images = 0

while np.sum(inertia[:num_of_images]) < objective:
    num_of_images += 1

print("")
print(f'Number of images to reach {objective*100}% of inertia: {num_of_images}')

Intertia:
[0.8206195  0.1065786  0.04214912 0.01383761 0.012708   0.00410717]

Number of images to reach 95.0% of inertia: 3

# Avec le bruit: 
img_b = load_image('Landsat_Kedougou/landsatKedougouC6bruit')
n_bands = 6

img = np.zeros((img_shape[0], img_shape[1], n_bands))
# Stacking up images into the numpy array:
for i in range(n_bands):
    if(i != 5):
        img[:, :, i] = load_image(path_img + str(i+1))
    else:
        img[:, :, i] = img_b

# Visualize them with imshow: 
n_rows = 2
plt.figure(figsize=(n_bands*10/n_rows, 20))

for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.imshow(img[:, :, i] / 255, cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

# Reshape the image for practicality: 
img = img.reshape(-1, n_bands)

# Standardization:
means = np.mean(img, axis=0)
stds = np.std(img, axis=0)

img = (img - means) / stds

# Compte the covariance matrix:
cov = np.cov(img, rowvar=False)

# Compute eigenvalues and vectors with eig:
eig_vals, eig_vecs = np.linalg.eig(cov)

# Lets sort the eigenvalues to be able to choose the components by importance:
idx = np.argsort(eig_vals)[::-1]
eig_vals = eig_vals[idx]
eig_vecs = eig_vecs[:, idx]

# Lets get the components in separate matrices, so we can display them easily:
img_principal = np.dot(img, eig_vecs)

plt.figure(figsize=(n_bands*10/n_rows, 20))

for i in range(n_bands):
    plt.subplot(n_rows,n_bands//n_rows,i+1)
    plt.imshow(img_principal[:,i].reshape(img_shape[0], img_shape[1]), cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

# Compute the percentage of inertia of each component:
inertia = eig_vals / np.sum(eig_vals)

print('Intertia:')
print(inertia)

objective = 0.95
num_of_images = 0

while np.sum(inertia[:num_of_images]) < objective:
    num_of_images += 1

print("")
print(f'Number of images to reach {objective*100}% of inertia: {num_of_images}')

Intertia:
[0.74345345 0.12236203 0.08137083 0.03492762 0.01373647 0.00414959]

Number of images to reach 95.0% of inertia: 4

def histogramme(image):
    hist = np.zeros(256, dtype=np.float32)
    shap = image.shape
    for j in range(shap[0]):
        for i in range(shap[1]):
            valeur = image[j,i]
            hist[int(valeur)] += 1
    return hist

path_img = 'SPOT/cam' 
n_bands = 3 
size = 512 

img = np.zeros((size, size, n_bands))

# Stacking up images into the numpy array:
for i in range(n_bands):
    img[:, :, i] = load_image(path_img + str(i+1))

n_rows = 1
plt.figure(figsize=(n_bands*10, 10))

# Use the "histogramme" function to have an intuition of pixel repartition on each image
for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.plot(histogramme(img[:,:,i]), color='black')

# Lieu du pixel caracteristique de l'eau
p1=[210, 150]
# Lieu du pixel caracteristique de la terre
p2=[50, 50]

from sklearn.cluster import KMeans

p1_vec = p1[0] * size + p1[1]
p2_vec = p2[0] * size + p2[1]

starting_points = np.array([[p1_vec], [p2_vec]])
img_vec = img.reshape((-1, n_bands))

kmeans = KMeans(n_clusters=2, init=starting_points, n_init=1)
labels = np.zeros((size, size, n_bands))

for i in range(n_bands):
  curr_img = img_vec[:, i].reshape(-1, 1) 
  kmeans.fit(curr_img)

  labels[:, :, i] = kmeans.labels_.reshape(size, size)

plt.figure(figsize=(n_bands*10, 12))
for i in range(n_bands):
    plt.subplot(1, n_bands, i+1)
    plt.imshow(labels[:,:,i], cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

path_img = 'Landsat_Kedougou/landsatkedougouC' 
n_bands = 6 
size = 256

img = np.zeros((size, size, n_bands))

# Stacking up images into the numpy array:
for i in range(n_bands):
    img[:, :, i] = load_image(path_img + str(i+1))


n_rows = 2
plt.figure(figsize=(n_bands*10, 35))

# Use the "histogramme" function to have an intuition of pixel repartition on each image
for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.plot(histogramme(img[:,:,i]), color='black')

# Fiding initial points for the water and land segmentation
plt.figure(figsize=(8, 8))
plt.imshow(img[:, :, 4], cmap='gray')

# Lieu du pixel caracteristique de l'eau
p1 = [180, 183]
# Lieu du pixel caracteristique de la terre
p2 = [172, 210]
# Lieu du pixel caracteristique de la terre 2
p3 = [185, 150]
# Lieu du pixel caracteristique de la terre 3
p4 = [220, 55]


plt.scatter(p1[0], p1[1], color='blue', s=10)
plt.scatter(p2[0], p2[1], color='red', s=10)
plt.scatter(p3[0], p3[1], color='purple', s=10)
plt.scatter(p4[0], p4[1], color='orange', s=10)

p1_vec = p1[0] * size + p1[1]
p2_vec = p2[0] * size + p2[1]
p3_vec = p3[0] * size + p3[1]
p4_vec = p4[0] * size + p4[1]

starting_points = np.array([[p1_vec], [p2_vec], [p3_vec], [p4_vec]])
img_vec = img.reshape((-1, n_bands))

kmeans = KMeans(n_clusters=4, init=starting_points, n_init=1)
labels = np.zeros((size, size, n_bands))

img_vec = img.reshape((-1, n_bands))

for i in range(n_bands):
  curr_img = img_vec[:, i].reshape(-1, 1) 
  kmeans.fit(curr_img)

  labels[:, :, i] = kmeans.labels_.reshape(size, size)

n_rows = 2
plt.figure(figsize=(n_bands*10//n_rows, 20))
for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.imshow(labels[:,:,i], cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

path_img = 'Landsat_Tarascon/landsattarasconC' 
n_bands = 8 
size = 512

img = np.zeros((size, size, n_bands))

# Stacking up images into the numpy array:
for i in range(n_bands):
    img[:, :, i] = load_image(path_img + str(i+1))


n_rows = 2
plt.figure(figsize=(n_bands*10, 35))

# Use the "histogramme" function to have an intuition of pixel repartition on each image
for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.plot(histogramme(img[:,:,i]), color='black')

kmeans = KMeans(n_clusters=3, n_init=10)
labels = np.zeros((size, size, n_bands))

img_vec = img.reshape((-1, n_bands))

for i in range(n_bands):
  curr_img = img_vec[:, i].reshape(-1, 1) 
  kmeans.fit(curr_img)

  labels[:, :, i] = kmeans.labels_.reshape(size, size)

n_rows = 2
plt.figure(figsize=(n_bands*10//n_rows, 20))
for i in range(n_bands):
    plt.subplot(n_rows, n_bands//n_rows, i+1)
    plt.imshow(labels[:,:,i], cmap='gray')
    plt.title('band '+str(i+1),fontsize=15)
    plt.axis('off')

SI 221-B

Bases de l’Apprentissage

TP "PCA & K-means"

MARTINS BRAGA Arthur

Données¶

PCA¶

K-means¶